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Abstract. Central schemes are frequently used for incompressible and compressible 
flow calculations. The present paper is the first in a forthcoming series where a new 
approach to a 2nd order accurate Finite Volume scheme operating on cartesian grids 
is discussed. Here we start with an adaptively refined cartesian primal grid in 3D 
and present a construction technique for the staggered dual grid based on L°°-Voronoi 
cells. The local rehnement constellation on the primal grid leads to a finite number 
of uniquely dehned local patterns on a primal cell. Assembling adjacent local patterns 
forms the dual grid. All local patterns can be analysed in advance. Later, running the 
numerical scheme on staggered grids, all necessary geometric information can instantly 
be retrieved from lookup-tables. The new scheme is compared to established ones in 
terms of algorithmical complexity and computational effort. 
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1 Introduction 


CFD simulations, especially in 3D, produce large amounts of data. In order 
to minimize memory requirements and computing time without sacrificing high 
spatial resolution in regions of interest modern numerical schemes are usually 
based on adaptive grids. Because of their simple structure and the ease of data 
access, cartesian grids are particularly popular. 

Staggered grids are a pair of meshes of the same computational domain whose 
nodes (in ID), edges (in 2D) and faces (in 3D) do not coincide. While the 
shape of a staggered grid is canonical on an uniform mesh it becomes rather 
complicated for an underlying adaptively refined grid, in particular in three 
space dimensions. In the present paper we suggest a construction technique for 
a staggered grid which is dual to an adaptive cartesian grid. 

Our focus of interest is the class of central schemes for hyperbolic conservation 
laws, which became quite popular during the last decade. In this case, staggered 
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grids circumvent the possibly costly numerical solution of Riemann problems. 
The prototype of all central schemes is the first-order Lax-Friedrich scheme 
[?]. Higher order extensions were first introduced by Nessyahu and Tadmor [?], 
and developed further by many others (see [?, ?, ?, ?, ?, ?] and the references 
therein). 

Let us briefly summarize our staggered grid construction. The cells of the dual 
grid are defined as the L°“-Voronoi regions around the vertices of the primal 
grid. Even though it is in principle simple to construct these Voronoi regions, it 
becomes a very complex task if the primal grid is adaptively refined, especially in 
three dimensions. Our approach, motivated by previous work on the 2D case [?], 
is based on a local decomposition of each primal cell into dual subcells. We call 
this decomposition of a primal cell the local pattern. The dual subcells are later 
assembled at the vertices. Since we use a cartesian grid with graded refinement, 
the number of different local decompositions is bounded. In [?], we used Polya 
theory to show that there are 227 combinatorially different local patterns in 3D. 
All these patterns can be analysed and stored in advance. Later, running the 
numerical scheme on staggered grids, all necessary geometric information can 
instantly be retrieved from lookup-tables. 

A popular alternative dual grid construction on adaptive cartesian grids is based 
on “diamond cells”, developed for example in [?, ?]. Although this geometrical 
description is very simple, we show in Section [5] that diamond grids increase the 
numerical cost of the finite volume scheme considerably. 

The paper is organized as follows: In Section [5] we investigate Voronoi decom¬ 
positions with respect to several norms and compare them with diamond dual 
grids. In order to explain the techniques as clearly as possible, we begin by con¬ 
structing the local patterns for 2D grids in Section[31 In Section^ we generalize 
the pattern construction to 3D grids. Section [5] presents first numerical applica¬ 
tions in 3D, discusses algorithmical complexity and poses some open questions 
to be treated in forthcoming papers. 


2 Problem setting and concept 


In this section we introduce the necessary notation, formulate properties of the 
primal and the corresponding dual grid and explain the basic idea of the dual 
grid construction. We discuss both the popular “diamond cell construction” as 
well as the Voronoi decomposition with respect to different norms in some detail 
and argue for our favorite choice. 

Throughout this paper we use the capital letters G, G, F, E and N for a 
grid, cell, face, edge resp. node, capital calligraphic letters to denote sets, like 
A/” = U V, subscripts p and d to distinguish between primal and dual objects, 
an asterisk * to denote objects from the corresponding staggered grid (the dual 
grid is staggered to the primal grid and the converse), subscripts i, j and k for 
node indices, and the superscript + to label local objects on single cells. 

Our primary interest of constructing dual grids is directed towards the numerical 
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solution of systems of conservation laws 

(1) Ut{x,t) + div F{u{x,t)) = 0 

in three space dimensions. Here, x G denotes the space variable, t the 
evolution time, u{x,t) G R"^ the solution vector in terms of m conservative 
variables, and F = the vector of the three directional flux functions 

/,;R-^R-,i = 1,2,3. 

The final goal is to implement a 3D-extension of the second order accurate 
central scheme proposed in [?]. A dimensionless formulation of the finite volume 
scheme reads as follows: 

Let V = v{x,t) be a cellwise smooth approximate solution of ([T]) with cell 
averages 

(2) j v{xX)dx. 

c 

The standard finite volume update is 

(3) J j F{v)-ndxdt 

t^dc 

Note that the flux integral should be replaced by a quadrature rule, and the 
approximate solution v needs to be extrapolated in time (for a fully discrete 
scheme). Since the solution v may be discontinuous at the cell boundary dC, 
one would have to replace the flux function F(v) by a Riemanii solver. This can 
be avoided when using staggered grids: now mean values on the timelevel 
are computed on cells C* of the corresponding staggered grid via 


(4) =t4 E 

' ' cnc*#0 


J v{x,F)dx+ J J F{v)-ndxdt 
cnc* ("cnac* 


h 


I 2 


As before, one needs quadrature rules for the integrals and time extrapolation 
for V. For the flux integral I 2 this is possible in such a way that F{v) needs to be 
computed only at points (x, t) where v is continuous [?, ?, ?]. Therefore, there is 
no need to solve Riemanii problems for central schemes. Still, the evaluation of 
the flux function may be the most expensive part of the finite volume scheme. 
Therefore, the numerical scheme should minimize the number of quadrature 
points in 12 - 

The cartesian primal grid Gp is supposed to be an adaptively refined 3D carte¬ 
sian mesh. We restrict ourselves to the regular subdivision of a hexahedron into 
eight child hexahedra. The grid adaptation is subject to a 1-level transition 
constraint between cells which share a common face or a common edge, see 
Figure [T] Any refinement technique which observes these rules would be ap¬ 
propriate, as for instance AMR [?], the cartesian refinement based on saturated 
error indicators [?, ?] or multi-scale analysis [?]. 
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Figure 1: 1-level transition constraint for the primal 3D-grid 


Given a primal grid Gp we are now looking for a corresponding staggered dual 
grid Gd- The numerical algorithms, and in particular those of finite volume 
type, require at least: 

• Gp and Gd are meshes for the same computational domain 0. 

• Interior nodes, edges and faces of Gp and Gd do not coincide. 

Moreover, the dual grid Gd should locally reflect the resolution of the primal 
grid Gp. With respect to the limitation of the timestep size by the CFL-number 
and the distance of discontinuities in the numerical solution one should also try 
to maximize the distance between faces of primal and dual cells. 

Before we present our new approach for the construction of the staggered grid, 
we first discuss the well known “diamond grids” (cf. [?, ?]). For reasons of clarity 
only the 2-dimensional figures are presented. All arguments are formulated also 
for 3-dimensional grids. The steps to construct dual cells read as follows: 

1. Construct edges of dual cells simply by connecting the midpoint of a primal 
cell G with all of C’s corners (see Figure [^(a)[ ). In 2D these edges subdivide 
G into four triangular dual cell parts, whereas in 3D the dual edges span 
12 dual faces which bound six pyramidal dual cell parts. 

2. Dual cells result by sticking together adjacent dual cell parts with a com¬ 
mon primal face (see Figure [2(b)[). 




(a) constructing dual edges (b) resulting dual diamond cells 


Figure 2: diamond-type dual cell construction in 2D 
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This construction technique is impressively simple, but suffers from excessive 
numerical cost: 

• Since every primal grid cell is subdivided into four (in 2D) resp. six (in 
3D) dual cell parts, and a dual cell is normally composed by two cell parts 
the dual grid contains approximately twice resp. thrice as much cells as 
the corresponding primal grid. 

• Eollowing the algorithm of Nessyahu and Tadmor from [?] two kinds of 
integration have to be performed in Q to determine the numerical solution 
on the next time level. The far more expensive part is I 2 , the integration of 
flux functions over the boundary of cells. Since the cellwise representation 
of the numerical data is in general discontinuous at the cell boundaries, 
the fluxes have to be evaluated in the interior of primal cells (e.g. in the 
center of gravity of the dual faces in the case of piecewise linear data. Its 
evaluation in the nodes of the dual cells is not possible since they often 
coincide with the boundary of primal cells). Since there are four (in 2D) 
resp. twelve (in 3D) dual faces per primal cell the number of necessary 
flux evaluations for a timestep onto the dual grid sums up to: 

^ fluxes(2D) = 4|Cp| 
fluxes(3D) = 12|Cp| 

• Due to the finer resolution of the dual grid the diameter of a dual cell is 
by factor smaller than the diameter of its corresponding primal cell. 
This leads to a restriction of the global timestep size in the explicit time 
integration procedure. 

Our less expensive alternative approach consists of the following principles: 

1. Construct dual cells always surrounding exactly one node Np G Af{Gp). 
Hence dual nodes Since in a cartesian grid there are approximately as 
many nodes as cells, primal and dual grid are of comparable size : 

# fluxes(2D) = 2|A^d| 2|Cp| 

# fluxes(3D) = 3|A/d| Ri 3|Cp| 

Moreover, the dual grid reflects to local resolution of the primal grid. 

2. The shape of the dual cells is determined by a Voronoi decomposition of 
the domain H respecting all nodes Np G Af{Gp). Nodes of dual cells now 
lie in the interior of primal cells. For the flux-integration the flux-functions 
can be evaluated there with respect to every spatial direction. Thus the 
numerical effort is essentially determined by the number of nodes in the 
dual grid (which is only slightly larger than the number of primal cells). 

The cartesian structure of the primal grid and its adaptation constraints (regular 
refinement, 1-level transition) lead to a rigidly prescribed location of the primal 
nodes. This allows the Voronoi decomposition to be performed locally on every 
cell Gp of the primal grid (as it will be outlined in the next chapter). Depending 
on the set of primal nodes on the boundary of a primal cell Gp, J\f'^{Gp) := 
J\f{Gp)r\Cp, we deduce Cp’s decomposition into local Voronoi regions. All these 
local regions together form a local pattern on Gp. Local patterns on adjacent cells 
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of the primal grid automatically match. All local Voronoi regions corresponding 
to a fixed node Np of the primal grid finally form the dual cell Cd{Np). 

Due to the primal grid structure and the 1-level transition constraint the number 
of essentially different local patterns is finite, and every possible decomposition 
can be analysed in advance. Later, the numerical scheme gets instant access to 
these predefined patterns and uses them in scaled and rotated copies. These 
copies are retrieved at run time and do not have to be stored. 


Voronoi regions: choice of the norm 

Given a set U C K" and a finite set of nodes Af = {Nk, fc = 1,..., m | Nk S U}, 
we decompose U into m Voronoi regions Vn^. by the following conditions: 

(i) uT=iyN. = ^ 

(ii) VNi n Vnj 

(hi) llx - Nk\\ < ||x - Nj\\,x e VN^yj k. 

Here denotes the interior of Vat^, and || • || an arbitrary norm on K". 

Let us now discuss the choice of an appropriate norm in 2D more in detail. The 
same arguments hold also for 3D. 

A Voronoi decomposition of U C refering to only two nodes Ni resp. Nj G U 
splits U into two Voronoi regions, ) (around Ni) and VAr.(Ar^) (around Nj), 

sharing a common separating polygon Stj. This separating polygon is the locus 
of the intersection points of circles with same diameter centered at Ni resp. Nj. 
The shape of the separating polygon depends both on the locus of Ni and Nj 
as well as on the used norm on In general, a Voronoi region Vn^ around the 
node Ni takes the form 

(5) Uvi = Pi Uvi(Ar„-) 


Respecting the ordinary || • || 2 -norm, any Sij is a straight line, hence Voronoi 
regions V/v, are always convex. In contrast, for the || • ||oo-norm, Sij consists, in 
general, of three straight lines, which are aligned with the cartesian axes (x-axis, 
y-axis) or the plane diagonals (of the xy-plane), see Figure 3(a) For special loci 
of Ni and Nj the separating polygon simplifies to a straight line (see Figure 
|3(b)[ ), or becomes not-unique in regions away of V and Nj (see Figure 3(c) I. 
For the latter case we choose the prolonged straight line from the unique part 
as the separating polygon. Moreover, Voronoi regions UjVi might happen to be 
non-convex. 


However, thanks to the cartesian structure of the primal grid the || • ||oo-iiorm 
even simplifies the construction of local Voronoi regions on cells of the primal 
grid. Compared to the corresponding construction for the Euclidean norm, 
boundary faces of Voronoi regions are now aligned only with the axes or the 
plane diagonals of the primal grid (cf. Figure [4(b)| vs. Figure [4(a)| ). We therefore 
favour the || • ||oo-norm. 
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Figure 3: Different configurations of two points, some || • ||oo-iiorm circles around 
these points, and the resulting separating polygons in 2D 



Figure 4: Comparison of Voronoi regions and separating polygons on a primal 
cell in 2D 


Local Voronoi construction 


Next, we inverstigate the locality of the Voronoi construction. With ® we get 
bounding boxes of local Voronoi regions Vvi on a primal grid cell Cp by 

(6) V^":=Cpn f| DCpCVv. 

NjeAfiCp) 


Af{Cp) denotes the set of corners of the cell Cp. Figure [5] illustrates the bounding 
boxes resulting by intersection (|6]) on a primal cell in 3D. Here, W is a corner 
node (fig. 5(a) I, a node in the midpoint of an edge (fig. |5(b)[ ) resp. a node in 
the midpoint of a face (fig. 5(c) I of Cp. 



Figure 5: Local Voronoi construction in 3D 


In order to get the final shape of a local Voronoi region Vn- on Cp it suffices to 
perform the intersection ([S]) for all nodes Nj G N'~^{Cp), i.e. respecting only the 
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nodes on the boundary of Cp. Shorter, we state 

Lemma 1. 

Cp n VNi = CpD Pi VMi(Nj) 

NjeM+iCp) 


Proof. Assuming lemma [T] were false. Then there is a node G Af\Af^{Cp) 
such that 0 U/Vi(Arfe) £ ^Ni which implies 

(7) 3xG with ||a; - A^fe||oo < \\x - A^ilU 

i.e. X should not be assigned VMi- 

Let S^Ni{.Cp) := ^||y-AfilU(l/)- Now, © implies 

( 8 ) N^G^nACp) 

With /S.X = diam(Cp, || • ||oo) we define ilAr^(Cp) := U Cp B^{y) and 




^Ni{Cp) \ Cp{Ni) if Ni is corner node of Cp 
^Ni{Cp) else 


where Cp{Ni) denotes that neighbour of Cp which shares only the node Ni, i.e. 

Cp n Cp(W) = N,. 


Suppose that we are in case 5(a) Let z G Cp{Ni) be fixed. Then we have 

’ O’ c 

\\y — Ni\\ao < \\y — ^||oo for any y G Honce there is no y G such that 
z G R||„_Ar,||(y), thus Cp{Ni) nf)Ni{Cp) = 0. 



the possibly occuring hanging nodes on Cp’s boundary. Hence 
(9) (Cp)nW\AA+(Cp) = 


□ 


i.e. there is no Nk fulfilling condition ([8]). This proves lemma [TJ 

Remark 1. The special treatment of the case 5(a)\ where Ni is a corner of Cp, 
is necessary since the 1-level grading condition over faces and edges still allows 
a 2-level transition over nodes. In this case hXNiiCp) and Af \ {Cp) would 
share a node of Cp{Ni), and ©) would not hold. 

Remark 2. The arguements in the proof do not depend on the spatial dimension 
of Cp. For a better understanding the reader might prefer to draw sketches of 
the proof in the 2-dimensional setting. 


3 Pattern construction in 2D 


In contrast to the local patterns from Figure 9(a) used for the dual grid con¬ 
struction in [?] we now want to create them by a clear geometrical rule. By 
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means of a local Voronoi decomposition respecting the maximum norm we as¬ 
sign the “volume” of a cell Cp of the primal grid the nodes Nk G on 

the boundary of Cp. These are the four corner nodes of Cp (they always exist) 
and possibly up to four nodes on the midpoint of Cp’s edges. These additional 
hanging nodes only exist if adjacent cells in the primal grid were refined. The 
local intersection of seperating polygons boils down to six different types of local 
Voronoi regions, cf. Figure [S) 



Figure 6: Resulting local Voronoi regions on a primal cell in 2D 


Obviously, all cuts of the cell Cp run parallel to the axes or the plane diagonals. 
We could therefore describe the steps for the construction of local Voronoi re¬ 
gions on Cp equivalently by the following algorithm: 

1. Subdivide the cell Cp into 16 congruent squares. 

2. Subdivide all squares containing the midpoint of Cp into two triangles. 
The cuts follow the face diagonals through Cp’s midpoint. On Cp we get 
8 triangles. 


Figure 7: Subdivision of a primal cell in 2D 

These squares and triangles, depicted in Figure [71 will be denoted as atoms 
of Cp. 

3. Finally, we assign these 20 atoms of Cp the nodes of J\f^{Cp) simply by 
calculating the || • ||oo-distance between the nodes and the atom’s center of 
gravity. All atoms being assigned the same node Nk G J\f'^(Cp) form the 
local Voronoi region VN^iCp) on Cp. The set of all local Voronoi regions 
on Cp forms the local pattern. Figure [8] shows all six essentially different 
local patterns in 2D. 



Figure 8: All essentially different local patterns in 2D 
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Assembling scaled and rotated copies of appropriate local patterns lead to the 


proposed in [?] there is now a unique correspondance between dual cells and 
primal nodes. This simplifies the adressing of data on the dual grid. 


dual grid depicted in Figure [9(b)[ Compared to the dual grid from Figure 9(a) 



(a) construction from [?| 



(b) Voronoi decomposition 


Figure 9: 2D cartesian primal grid (solid) and corresponding dual grid (dashed) 


4 Pattern construction in 3D 


Now we extend the steps to create local patterns described in the previous sec¬ 
tion to 3-dimensional grids. Applying the concept of constructing local Voronoi 
regions in 3D leads to an only slightly more complicated description than in 
2D. Again, we assign the volume of a primal cell Cp the nodes on its boundary, 
i.e. all Nk S J\f'^{Cp). These are the eight corner nodes of Cp, and possibly 
up to 18 nodes on the midpoint of Cp’s six faces and 12 edges. These hanging 
nodes occur if an adjacent cell, sharing a common face resp. edge with Cp, was 
refined. For a general location of two nodes Ni and Nj in the separating 
surface between the Voronoi regions Vv^(Ar^) and consists of up to seven 

planar patches, whose normals are aligned with the axes {x-, y-, z-axis) or the 
plane diagonals {xy-, yz-, za:-plane). 

Fortunately, the location of nodes in the primal grid restricts the variety of 
occuring separating surfaces to the six cases of Figure [TOl Finally, on a primal 
cell only scaled, translated and rotated versions of these cuts have to be executed 
to get the local Voronoi regions. 



Figure 10: Occuring separating surfaces in 3D 


Hence the whole construction can again be described as a subdivision algorithm: 


1. Subdivide the hexahedral cell Cp into 64 congruent cubes. 
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(a) subdivision of a pri¬ 
mal cell 



(b) local Voronoi regions 



(c) non-convex Voronoi 
region 


Figure 11: Construction of local Voronoi regions in 3D 


2. Subdivide all cubes containing a midpoint of Cp’s faces into two prisms. 
The cuts follow the face diagonals through the face midpoints. On Cp we 
get 8 X 6 = 48 prisms. 

3. Subdivide all cubes containing the central point of Cp into six tetrahedra. 
The corresponding three cuts equal the diagonal cuts on their adjacent 
cubes from step [H This results in 6 x 8 = 48 tetrahedra. 


These cubes, prisms and tetrahedra will again be denoted as atoms of Cp. Figure 
|ll(a)| shows a cell Cp and some of its atoms. Next we assign these 64-1-48 — 24-1- 
48 —8 = 128 atoms of Cp the nodes oiAf'^{Cp) by calculating the || • ||oo-distance 
between the nodes and the atom’s center of gravity. All atoms being assigned 
the same node Nk S Af~^{Cp) form the local Voronoi region VN^iCp) on Cp. 
Figure 11(b) shows a cell Cp, the set {Cp)\N{Cp) (i.e. the corners of Cp left 

out) and a few corresponding local Voronoi regions. Let us mention that local 
Voronoi regions might happen to be slightly non-convex (region in the center of 
Figure [Tl(c)[ ). But this should not effect the feasibility of our approach, and we 
do not expect additional problems in constructing a second order central finite 
volume scheme by the occurrence of these dual cells. 

The set of all local Voronoi regions on Cp form the local pattern. The shape 
of these three-dimensional local patterns on the boundary faces of a primal cell 
matches the two-dimensional patterns from Figure |51 


Counting local patterns 

Since the number of different distributions of nodes on the boundary of a pri¬ 
mal cell is finite, all local patterns can be assembled and stored in advance. 
Later, running the numerical scheme, these patterns are used in scaled and 
rotated copies. To determine the number of essentially different local pat¬ 
terns (that do not arise from each other by rotation or reflection) we used 
tools of the Polya theory (see [?], [?]). All information for the counting prob¬ 
lem is provided by the group of symmetries 'D{W) of the cube 211. Analysing 
how the elements of 'D{W) act on the set of QR’s edges and faces leads to 
the number of 227 essentially different local patterns. For details we refer to 
[?]. The list of all 227 patterns (some of them with figures) can be found on 
http://www.igpm.rwth-aachen.de/wolfrain/local_patterns.html 
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They can all be assembled, analyzed and stored in a lookup-table. Every local 
refinement situation can be identified with exactly one reference pattern. This 
reference pattern as well as an appropriate mapping (i.e. permutation matrix) 
onto the really occuring pattern are a priori known. 


5 Applications 


At the current stage, our implementation of the finite volume scheme includes 
the primal and dual grid generation, the integration of piecewise linear func¬ 
tions over volumes and boundary faces of cells, but still lacks some essential 
components (like data reconstruction and limiting, treatment of divers bound¬ 
ary conditions) to end up with a second order accurate non-oscillatory method. 
Nevertheless, we have verified our dual-grid approach with standard test cases. 

After summarizing the algorithmical demands and presenting our integration 
strategy, we address the data access, describe our test cases and finally compare 
the expected numerical effort of our method with the diamond grid approach 
and the established non-staggered HLL-scheme. 


Algorithmical demands 

The evaluation steps in o are twofold. Part Ii calls for an integration of 
the piecewise linear function u" over a bounded volume and can be evaluated 
exactly. Part I 2 is a bit more intricate. For an illustration we refer to Figure 
M The normal part of a non-linear flux-function F has to be integrated in time 
over an interval At = — F and in space over a polygonally bounded planar 

cell face p. In general, a face p of a cell C* consists of several polygonal parts 
pi in adjacent cells Ci of the corresponding staggered grid. 

Any second order accurate quadrature formula for I 2 should evaluate F in time 
and space where its argument v = v{x,t) is uniquely defined. For the temporal 
part of I 2 we use the midpoint rule. For that reason v has to be extrapolated to 
approximate As time evolves, characteristics spread out from the 

boundary into the interior of each cell. In order to avoid the solution of local 
Riemann problems at the cell interfaces the spatial integration in I 2 should 
choose quadrature points away from the boundaries which are, at time 
not yet reached by the characteristics. The distance of these points to the cell 
boundary limit the size of the global timestep At = — t". 

Though admissible, the choice of quadrature points in the center of gravity of 
every polygonal part pi is disadvantageous: 

• In a cartesian 3D-grid there about three times more faces than cells. A 
polygonal face contains of mostly four polygonal parts. We would have to 
evaluate one flux function on every polygonal parts. 

• The distance of the center of gravity of a polygonal part pi to the boundary 
of the cells Ci of the corresponding staggered grid is rather short. This 
limits the global timestep size unnecessarily. 
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In our algorithm we implemented a generalized trapezoidal rule on polygons. 
The flux function F is evaluated at the nodes of a polygonal face. Later, 
these values are scaled with barycentric weights respecting the face geometry to 
achieve a second order accurate integration rule. This approach promises to be 
more competitive: 

• In a cartesian 3D-grid the number of nodes and the number of cells is 
similar. We evaluate all three directional flux functions at every node 
Di. The numerical effort is by factor si 4 less than the method described 
above. 

• Since the distance of nodes of a polygonal face p of a cell C* to the 
boundary of cells Ci from the corresponding staggered grid is larger than 
the distance from the pi’s center of gravity to the boundary the global 
timestep size is less limited compared to the method above. 



Figure 12: Compound boundary faces of dual cells 


Data access 

For fast access all geometric and numerical data, both for the primal and the 
dual grid, are stored in (some of them level-linked) hashmaps. On every cell C 
of the grid G the finite volume scheme needs to know 

1. the volume of C 

2. its neighbouring cells on the grid G (for a higher order reconstruction) 

3. the common faces of C and its neighbours, including their normals and 
their area, as well as the nodes on the grid G that form these faces (to 
determine fluxes over faces by an appropriate quadrature rule) 

On the cells of the cartesian primal grid Gp such information is quite easily ac¬ 
cessible by index shift operations. In contrast, for a cell Gd of the dual grid Gd 
this knowledge is scattered over all local patterns contributing to the construc¬ 
tion of Gd and has to be “collected” during a complete traversal of the primal 
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grid Gp. Hence finding the appropriate local pattern on a primal cell is an often 
used operation and has to be fast. 

The local pattern on a primal grid cell Gp is determined by the occurrence of 
hanging nodes in the midpoints of Gp’s 12 edges and six faces (cf. Section S) 
These information form an 18-bit key which allows instant access to the corre¬ 
sponding reference pattern (one of 227, cf. Section 0]) as well as the appropriate 
permutation matrix from a predefined hashmap. This hashmap can easily be 
initialized in a pre-roll step by applying each of the 48 self-mappings of the cube 
to all of the 227 essentially different local refinement constellations. Finally, one 
gets 6210 hashmap entries. 

A local pattern consists of a certain number of local Voronoi regions , each 
of them composed of several atoms (cubes, prisms and tetrahedra, see Figure 


1. the volume of Vn^. equals the sum of the volumes of its atoms 

2. neighbouring Voronoi regions on a pattern refer to neighbouring cells on 
the dual grid Gd 

3. common faces of neighbouring Voronoi regions are the union of common 
faces of their atoms. Thus, the normals, areas and the face-forming nodes 
are known by construction of the local pattern. 

Whereas on the primal grid Gp all cells Gp are cubes and common faces to 
neighbouring cells are always square-shaped, there is a wider range of different 
cell- and face-types on the dual grid Gd, see Figure [T3( a) | for a rough impression. 


11(a) I. We get all information listed above locally: 


Example I: Volume integration 

In order to test the presented dual grid construction and to roughly estimate the 
numerical cost of a staggered grid scheme, a volume integration over a subset 
of the dual grid Gd was performed. The primal grid has been constructed by 
adaptively resolving the isosurface of a tilted elliptical paraboloid inside the unit 
cube, see Figure [r3(b)[ 

We verified our algorithm by means of the Gauss integral theorem 

/ div X) dx = — X) -n dx 

Jh JaCi 

by choosing d(x) = ^x, Vx S a subdomain U C fl, and integrating over all 
cells Gd C U of the dual grid Gd 

1. by summing up the volumes of the involved Voronoi regions 

2. by integrating i(x, n/) over all common faces / of neighbouring Voronoi 
regions (in view of flux integration over cell boundaries). Here n/ denotes 
the oriented normal of /. 

We performed the volume integration described above for several resolutions of 
the primal grid Gp and examined 
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(a) Close-up view on the dual grid. Lines 
represent the primal grid. 


(b) primal test-grid with resolved isosur¬ 
face on level 7, 112470 cells 


Figure 13: Primal and dual grid 


• the number of different occuring local patterns: 
it leveled off at about 40 

• the frequency of occurrence of these patterns: 

approximately 80% of the cells Cp on the primal grid refer to the most sim¬ 
ple pattern (no additional node on the boundary of Cp), whose numerical 
cost are low 

• the expected number of flux evaluations for the staggered grid solver by 
counting the number of nodes in the dual grid 

To estimate the profit of our approach, the numerical effort was compared with a 
cheap iioii-staggered finite volume scheme where flux integrals are approximated 
by the Harten-Lax-van Leer (HLL) Riemann solver (cf. [?]). Results are listed 
in table [T] below. The numerical cost of both schemes are similar. 


level 

primal cells 

dual cells 

fluxes non-staggered 

fluxes staggered 

4 

568 

1019 

3810 

4779 

5 

4502 

6969 

30126 

33717 

6 

24781 

34710 

163899 

161001 

7 

112470 

149011 

737283 

690420 

8 

493368 

630364 

3210945 

2906742 

9 

2156547 

2684250 

13944345 

12220641 


Table 1: Comparison between staggered and non-staggered approach 


Example II: Rotating Cone 

Not yet having incorporated linear data reconstruction and limiting in our 
method (what would be essential for a second order non-oscillatory scheme) 
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we run a standard linear advection problem, the “Rotating cone”. The compar¬ 
ison to the diamond-grid approach as well as to the non-staggered HLL-scheme 
in terms of numerical effort confirms our topological statements from Section [5] 
and encourages for the further work. 


Initial data & numerical solution 


As initial data for the rotating cone advection problem we resolve a smooth 
function / with compact support over an octant of the surface dB of a ball 
B = B{c,R) (cf. Figure 14(a)) The ball B is centered at c = {cx,Cy,Cz) = 
(0.6, 0.3,0.2) and has a radius R= 1/4. The function / is defined as follows: 


= (x- Cxf + (y - CyY + {z- Czf 

q = 4|l-r^| 

{ 1 — 2q‘^ for q < 1/2, x,y,z >0 
2(g-l)2 for 1/2 < g < l,a;,?/,z > 0 
0 else 


This cone is now rotated around c in the plane spanned hy v = (1,0,0)* and 
w = (0,1,0.5)*. The advection performs up to time t = 7r/4 at constant angle 
velocity 1. The exact solution is depicted in Figure [14(b) [ Since our scheme 
is still only first order accurate the numerical solution suffers from immoderate 
smearing (cf. Figure 14(c)). 


Grid specifications &: numerical effort 

Nevertheless, the resulting grids offer an opportunity to estimate and to com¬ 
pare the numerical effort of our scheme. The method to count the number of 
necessary flux evaluations in two successive timesteps is outlined in table [5] and 
performed for the grids from Figure [T4(c)| in table [3] 


timestep 

Diamond 

HLL 

Voronoi 

primal —>■ dual 
dual —>■ primal 

12 ^primal cells 

1 ^/Iprimal faces 

2 ^primal faces 

2 ^primal faces 

3 #dual nodes 

3 ^primal nodes 


Table 2: Counting flux evaluations in two successive timesteps 


//: Flux evaluations 


Gp —^ Gd 

Gd —>■ Gp 

E#/ 

Diamond 

1.35 M 

0.35 M 

1.69 M 

HLL 

0.69 M 

0.69 M 

1.38 M 

Voronoi 

0.46 M 

0.36 M 

0.82 M 


Grid specifications 

primal cells 

112106 

primal faces 

345564 

primal nodes 

121561 

dual nodes 

152630 


Table 3: Grid specifications and numerical effort 

Obviously, in an adaptively octtree-refined cartesian 3D-grid the number of 
cells on the (locally) deepest refinement level, i.e. those without hanging nodes 
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(a) Initial data at time 0 



(b) Exact solution at time 71/4 



(c) Numerical solution at time 7r/4. 


Figure 14: Rotating cone with first order scheme 


(we call it property %), outranges the number of coarser cells by far. In gen¬ 
eral, discontinuities in the numerical solution are resolved by (possibly thin) 
2-dimensional manifolds. For a rough argument, let be the number of grid 
cells and Ar the number of cells with property "H on a grid with at most L re¬ 
finement levels. Resolving a discontinuity by refining k cells on refinement level 
L effects the numbers A and N as follows: each of the k cells generates 8 child 
cells, the original cell itself disappears, thus N^+i « + 7k. The refinement 

of a cell inserts hanging nodes on its coarser neighbour cells. Since all k cells 
are supposed to form a thin surface, each of them “destroys” the property TL 
only in the two face-neighbours in the surface’s normal direction. The refined 
cell itself contributed most likely to A^. Every child cell on level L -I- 1 obvi¬ 
ously has property H. Thus Al+i ~ Al -I- 5k. For k ^ 1 one gets the ratio 
r = A/N Ri 5/7. If the k cells on the fine grid level rather form a volume than 
a thin surface, the ratio r gets even closer to 1. 

In this regard, the ratio of 94% in our example from Figure |14(c)| does not 
surprise. For all these cells the according local pattern is the easiest possible, 
only containing one dual node which leads to only half as much local costs than 
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incurred by the compared HLL-solver, and only one fourth in comparison to the 
diamond cell approach. 


6 Conclusions & extensions 

In the current paper basic steps of the implementation of Nessyahu and Tad- 
mor’s staggered scheme [?] on adaptively refined cartesian grids in 3D were 
introduced. The construction of dual grid cells as Voronoi regions respecting 
the maximum-norm simplifies their geometrical description and, in the end, the 
implementation. Substantial combinatorial investigation in advance allows our 
scheme in return to exploit local geometrical structure of the grid at run-time 
to save a lot of numerical effort. 

Though up to now only first order, the extension of the staggered grid ap¬ 
proach to a higher order 3D finite volume scheme is in process. It concerns the 
piecewise linear reconstruction and limitation of cell- and flux-values to end up 
with a second order scheme, and the data handling for large 3D simulations. 
Still challenging are the incorporation of obstacles in the computational do¬ 
main and the treatment of divers boundary conditions, especially when general 
CAD-geometries (and not only cartesian grid cells) need to be incorporated. 
Additionally desirable would be the parallelization of the staggered grid scheme 
as well as the local adaptivity in time. 



